#include "ProcessFunctions.h"

bool array2txt(const std::string &H, const char *cFilename) // 针对字符数组的输出函数
{
	std::ofstream outfile;
	outfile.open(cFilename, std::ios::out);
	if (!outfile.is_open())
	{
		std::cout << "Open file failure\n";
		outfile.close();
		return false;
	}

	for (int i = 0; H[i] != '\0'; i++) //遍历数组写入txt
	{
		outfile << H[i];
	}
	outfile.close();
	return true;
}

float *outHistogramData(HXLBMPFILE &src, const char *cFilename) //将绘制“灰阶-概率”直方图的数据保存到txt文件中
{
	float H[256] = {0}; // 8bit的灰度图，256灰阶
	int i, j;
	float *result = (float *)malloc(256 * sizeof(float));
	int total = src.imageh * src.imagew;
	for (i = 0; i < src.imageh; i++)
		for (j = 0; j < src.imagew; j++)
		{
			H[src.pDataAt(i)[j]] += 1.0; //对每个灰阶进行计数
		}
	for (i = 0; i < 256; i++)
	{
		H[i] /= (float)total; //计算图像概率p(Sk)
		result[i] = H[i];
	}

	if (!array2txt(H, cFilename))
	{
		printf("outHistogramData error!\n");
		exit(0);
	}
	return result;
}

float *HistEqualize(HXLBMPFILE &src, HXLBMPFILE &des, const char *cFilename) // 将src图像的直方图均衡化，并输出结果图像des
{
	float H[256], count_origin[256] = {0}, T[256] = {0}, count_new[256] = {0};
	// count_origin用来记录原图像各灰阶个数，count_new是均衡化之后的结果
	float *result = (float *)malloc(256 * sizeof(float)); // 返回均衡化后的灰阶-概率矩阵
	memset(H, 0, 256 * sizeof(float));
	int i, j;
	int total = src.imageh * src.imagew;

	for (i = 0; i < src.imageh; i++)
		for (j = 0; j < src.imagew; j++)
		{
			H[src.pDataAt(i)[j]] += 1.0;			// 先对灰阶计数
			count_origin[src.pDataAt(i)[j]] += 1.0; // 保留原图的计数结果
		}

	for (i = 0; i < 256; i++)
		H[i] /= (float)total; // 计算原图像概率p(Sk)

	for (i = 0; i < 256; i++)
		for (j = 0; j <= i; j++)
		{
			T[i] += H[j];	  // 求得均衡化之后的tk
			result[i] = T[i]; // 如果需要作图的那么result=count_new，如果需要均衡化result=Tk
		}
	if (!array2txt(T, "./output/tk_data.txt"))
	{
		printf("outHistogramData error!\n");
		exit(0);
	}

	for (i = 0; i < 256; i++)
	{
		T[i] = int((256 - 1) * T[i] + 0.5); // 用式tk=int((L-1)*tk+0.5)将tk扩展至[0,L-1]范围，得到灰度级
	}

	for (i = 0; i < src.imageh; i++)
		for (j = 0; j < src.imagew; j++)
		{
			des.pDataAt(i)[j] = T[src.pDataAt(i)[j]]; // 将新的值赋予给输出图像
		}

	for (i = 0; i < 256; i++)
	{
		count_new[int(T[i])] += count_origin[i]; // 根据新的灰度级求新的计数
	}

	for (i = 0; i < 256; i++)
	{
		count_new[i] /= (float)total; // 计算新图像概率p(tk)
									  // result[i] = count_new[i];
									  // 如果需要作图的那么result=count_new，如果需要均衡化result=T
	}

	if (!array2txt(count_new, cFilename))
	{
		printf("outHistogramData error!\n");
		exit(0);
	}
	des.SaveBMPFILE("./output/QI.bmp");
	return result;
}

float *HistEqualize(HXLBMPFILE &src, HXLBMPFILE &des, float *hist_data, const char *cFilename) //已知灰阶-概率矩阵，做直方图均衡化
{
	float count_origin[256] = {0}, T[256] = {0}, count_new[256] = {0};
	// count_origin用来记录原图像各灰阶个数，count_new是均衡化之后的结果
	float *result = (float *)malloc(256 * sizeof(float)); // 返回均衡化后的灰阶-概率矩阵
	int i, j;
	int total = src.imageh * src.imagew;

	for (i = 0; i < 256; i++)
	{
		count_origin[i] = hist_data[i] * total; // 求原图中的灰阶计数
	}

	for (i = 0; i < 256; i++)
		for (j = 0; j <= i; j++)
		{
			T[i] += hist_data[j]; // 求得均衡化之后的tk
			result[i] = T[i];	  // 如果需要作图的那么result=count_new，如果需要均衡化result=T
		}
	if (!array2txt(T, "./output/tk_data.txt"))
	{
		printf("outHistogramData error!\n");
		exit(0);
	}

	for (i = 0; i < 256; i++)
	{
		T[i] = int((256 - 1) * T[i] + 0.5); // 用式tk=int((L-1)*tk+0.5)将tk扩展至[0,L-1]范围，得到灰度级
	}

	for (i = 0; i < src.imageh; i++)
		for (j = 0; j < src.imagew; j++)
		{
			des.pDataAt(i)[j] = T[src.pDataAt(i)[j]]; // 将新的值赋予给输出图像
		}

	for (i = 0; i < 256; i++)
	{
		count_new[int(T[i])] += count_origin[i]; // 根据新的灰度级求新的计数
	}

	for (i = 0; i < 256; i++)
	{
		count_new[i] /= (float)total; // 计算新图像概率p(tk)
									  // result[i] = count_new[i];		//如果需要作图的那么result=count_new，如果需要均衡化result=T
	}

	if (!array2txt(count_new, cFilename))
	{
		printf("outHistogramData error!\n");
		exit(0);
	}
	des.SaveBMPFILE("./output/QI.bmp");
	return result;
}

float calculateA() // 计算 y=a(1/4-(x-1/2)^2)中的a
{				   // 一种比较直观且简单的算法是分割成多份然后求面积和，最后用1/sum即可求得a（因为概率和应为1）
	int i;
	float x[1000]; // 定义域[0,1]分割成0.001为步长的1000份
	float y[1000]; // 值域
	float sum = 0;
	for (i = 0; i < 1000; i++) // 由于c语言只能用整数下标，因此只能将x取值存为数组
	{
		x[i] = i * 0.001;
	}
	for (i = 0; i < 1000; i++) // 计算x对应的y值
	{
		y[i] = (0.25 - pow((x[i] - 0.5), 2)); // y=(1/4-(x-1/2)^2)
	}
	for (int i = 0; i < 1000; i++) // 计算面积和
	{
		sum += y[i] * 0.001;
	}
	return (1 / sum);
}

float *getTargetHist(const char *cFilename) // 由函数y=a(1/4-(x-1/2)^2)得到目标直方图，并返回灰阶-概率数组
{
	float a = calculateA(); // 先计算得到系数a
	float x[256], y[256];
	float *result = (float *)malloc(256 * sizeof(float));
	int i;
	for (i = 0; i < 256; i++)
	{
		x[i] = i * pow(2, -8); // 我们需要将函数扩展至0-255的范围，因此先划分
	}
	for (i = 0; i < 256; i++)
	{
		y[i] = a * (0.25 - pow((x[i] - 0.5), 2));
		y[i] /= 256;
		result[i] = y[i];
	}
	if (!array2txt(y, cFilename))
	{
		printf("outHistogramData error!\n");
		exit(0);
	}
	return result;
}

int *HistSpecificate(HXLBMPFILE &src, HXLBMPFILE &des, float *hist_data, const char *cFilename) //直方图规定化
{
	int i, j, index = 0;
	float *SQH_tk = HistEqualize(src, des, "./output/SQH.txt");			   //获取原图像的均衡化tk数组
	float *EQH_tk = HistEqualize(src, des, hist_data, "./output/EQH.txt"); //获取目标图像的均衡化tk数组
	int T[256];															   //灰阶对应关系矩阵
	int *result = (int *)malloc(256 * sizeof(int));
	float min, temp;
	for (i = 0; i < 256; i++)
	{
		min = fabs(SQH_tk[i] - EQH_tk[0]); // 假设最小值
		for (j = index; j < 256; j++)	   // 由于两个数组都是递增的，所以我们每次从上一次的结果继续即可
		{
			temp = fabs(SQH_tk[i] - EQH_tk[j]); // 依次做差找到最小值
			if (temp < min)
			{
				min = temp;
				index = j; // 找到最小值的下标
			}
			if (SQH_tk[i] < EQH_tk[j])
				break;
		}
		T[i] = index; // 保存映射结果
		result[i] = index;
	}
	if (!array2txt(T, "./output/HistSpecificate.txt"))
	{
		printf("outHistogramData error!\n");
		exit(0);
	}

	for (i = 0; i < src.imageh; i++)
		for (j = 0; j < src.imagew; j++)
		{
			des.pDataAt(i)[j] = T[src.pDataAt(i)[j]]; // 将新的值赋予给输出图像
		}
	des.SaveBMPFILE("./output/DI.bmp");
	outHistogramData(des, cFilename);
	return result;
}